evapol Function

      pure real(RKIND) function evapol(tu,nu,tv,nv,c,rad,x,y) result(e_res)

      !  calling sequence:
      !     f = evapol(tu,nu,tv,nv,c,rad,x,y)
      !
      !  input parameters:
      !   tu    : real array, length nu, which contains the position of the knots in the u-direction.
      !   nu    : integer, giving the total number of knots in the u-direction
      !   tv    : real array, length nv, which contains the position of the knots in the v-direction.
      !   nv    : integer, giving the total number of knots in the v-direction
      !   c     : real array, length (nu-4)*(nv-4), which contains the b-spline coefficients.
      !   rad   : real function subprogram, defining the boundary of the approximation domain. must be
      !           declared external in the calling (sub)-program
      !   x,y   : the co-ordinates of the point where f(x,y) must be evaluated.
      !
      !  output parameter:
      !   e_res : the value of f(x,y)
      !
      !  other subroutines required:
      !    bispev,fpbisp,fpbspl
      !
      !  references :
      !    de boor c : on calculating with b-splines, j. approximation theory 6 (1972) 50-62.
      !    cox m.g.  : the numerical evaluation of b-splines, j. inst. maths applics 10 (1972) 134-149.
      !    dierckx p. : curve and surface fitting with splines, monographs on numerical analysis, oxford
      !                 university press, 1993.
      !
      !  author :
      !    p.dierckx
      !    dept. computer science, k.u.leuven
      !    celestijnenlaan 200a, b-3001 heverlee, belgium.
      !    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
      !
      !  ..scalar arguments..
      integer,     intent(in) :: nu,nv
      real(RKIND), intent(in) :: x,y
      !  ..array arguments..
      real(RKIND), intent(in) :: tu(nu),tv(nv),c((nu-4)*(nv-4))
      !  ..user specified function
      procedure(boundary) :: rad

      !  ..local scalars..
      integer :: ier
      integer, parameter :: liwrk = 2, lwrk = 8
      real(RKIND) :: u(1),v(1),r,f(1),dist
      !  ..local arrays
      real(RKIND) :: wrk(lwrk)
      integer :: iwrk(liwrk)
      !  ..
      !  calculate the (u,v)-coordinates of the given point.
      u    = zero
      v    = zero
      dist = x**2+y**2
      if (dist>zero) then
         v(1) = atan2(y,x)
         r    = rad(v(1))
         if (r>zero) &
         u(1) = min(sqrt(dist)/r,one)
      endif

      ! evaluate s(u,v)
      call bispev(tu,nu,tv,nv,c,3,3,u,1,v,1,f,wrk,lwrk,iwrk,liwrk,ier)

      ! Return scalar result
      e_res = f(1)
      return
      end function evapol